Time-dependent versus static quantum transport simulations beyond linear response 
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To explore whether the density-functional theory non-equilibrium Green's function formalism 
(DFT-NEGF) provides a rigorous framework for quantum transport, we carried out time-dependent 
density functional theory (TDDFT) calculations of the transient current through two realistic molec- 
ular devices, a carbon chain and a benzenediol molecule inbetween two aluminum electrodes. The 
TDDFT simulations for the steady state current exactly reproduce the results of fully self-consistent 
DFT-NEGF calculations even beyond linear response. In contrast, sizable differences are found with 
respect to an equilibrium, non-self-consistent treatment which are related here to differences in the 
Kohn-Sham and fully interacting susceptibility of the device region. Moreover, earlier analytical 
conjectures on the equivalence of static and time-dependent approaches in the low bias regime are 
confirmed with high numerical precision. 



A first principles description of quantum transport at 
the atomic scale is usually based on the non-equilibrium 
Green's function (NEGF) technique [1], which reduces to 
the traditional Landauer transmission formalism [2 for 
coherent transport. The device Green's function is con- 
structed in an ad-hoc manner from the Kohn-Sham (KS) 
single-particle Hamiltonian of ground state Density Func- 
tional Theory (DFT) including the effects of the leads 
through suitable self-energies [3 . While such methods 
provide qualitative results to understand the transport 
behavior, there are cases where the theoretical steady 
state current differs from experimental results by orders 
of magnitude [4 . One reason for the discrepancy is that 
the evaluated transmission function exhibits resonances 
at the KS single particle energies, which are in general 
not physical. 

Recently, several approaches based on time- dependent 
DFT (TDDFT) have been put forward [HE], not only 
to treat genuine dynamical problems like AC transport 
or transient effects, but also to remedy the problematic 
situation in the steady-state case. For optical absorption 
spectra, TDDFT excited state energies improve signifi- 
cantly on single-particle energy differences that are ob- 
tained from the static ground state KS equations. Since 
the screening properties of the device depend strongly 
on the optical gap, one would expect also an improved 
description of the current. 

In fact, dynamical corrections to the conventional 
DFT-NEGF picture have been predicted in the regime of 
linear response 0. Interestingly, these corrections were 
shown to vanish for approximate xc functionals which are 
local both in time and space. Such adiabatic functionals, 
like the adiabatic local density approximation (ALDA 
a.k.a. TDLDA), are quite common in many implementa- 
tions of TDDFT for optical spectra due to their numerical 
simplicity. As shown by Stefanucci et. al. [8 and Zheng 
et. al. [6 , TDLDA leads to the well-known Landauer- 
Biittiker formula for the steady state current if it can 
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FIG. 1. (a) Schematic device setup, (b) Seven- membered 
carbon chain between Al nanowires in (001) direction. Fron- 
tier carbon atoms in hollow position of the Al surface at 1.0 
A distance. Bond lengths (A): C-C = 1.32, Al-Al = 2.86. (c) 
1,4-benezenediol molecule inbetween Al nanowires. 

finally be reached. It is thus argued that TDDFT and 
DFT-NEGF should result in the exact same steady state 
current if the adiabatic approximation is adopted for the 
exchange-correlation functional. We thus encounter a 
paradox. 

A related open question is whether sizable differences 
between TDDFT and DFT-NEGF may occur for local 
functionals beyond linear response. The interest is here 
in the full I-V characteristics of the device including the 
bias region around molecular resonances, which is often 
used as fingerprint of the contacted molecule in experi- 
ments. This regime is difficult to access with analytical 
methods but amenable to numerical simulations. Beyond 
linear response, more than one solution for steady state 
current is conceivable, and this was observed for model 
systems [9 . In such a case, the transient current may set- 
tle into different steady currents depending on how the 
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bias voltage is turned on. Also, it was argued that per- 
sistent oscillating currents can exist and a steady state 
may never be reached. This is because a device may have 
bound states with vanishing line widths and correspond- 
ing infinite life times. The implications of such bound 
states for open system were investigated for a model sys- 
tem [To], and the current through the model device was 
found to oscillate for as long as the simulation time. To 
address these issues in this letter, we use a recently devel- 
oped TDDFT method for open systems [6 and compare 
results with those from DFT-NEGF calculations. 

Methodology. — For a two-terminal setup, the entire 
system consists of a device region (D), and semi-infinite 
left (L) and right (R) electrodes as shown in Fig [l^. The 
region where electron-scattering events take place, is 
thus an open electronic system. 

To obtain the reduced electronic dynamics of the de- 
vice, we focus on the equation of motion (EOM) of tro, 
the reduced single-electron density matrix for D, 




(1) 



FIG. 2. TDLDA transient current of a) carbon chain and b) 
benzenediol for exponential turn-on of the bias voltage with 
different time constants T and Vb = ^V. 



where /id(^) denotes the KS Fock matrix of D, and the 
square bracket on the right-hand side denotes a commu- 
tator. Qa is the dissipation term due to the ath electrode 
and based on the Keldysh NEGF formalism [1 , it can be 
expressed as. 



Qa{t) 



dr 



G&(t,r)S<(r,t) + G<(t,r)S«(r,t) 



H.c. 
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where and denote the time domain retarded and 
lesser Green's function in the device region, and X)^ and 
T,^ stand for the advanced and lesser self-energies of lead 
a. 

The transient electric current through the interface Sa 
(the cross section separating region D from the ath elec- 
trode) can be evaluated via [6 : 



Ia{t) = -tT [Qa{t)]- 



(3) 



If the bias potentials approach a constant value asymp- 
totically {V^ = Va{t oo)), a steady state may de- 
velop. It has been shown that the steady state current 
can then be recast in the form O [8] : 

T{e)=tv[G^^{e)Tn{e)GUe)TL{e)], (4) 

where fa are Fermi distribution functions and Ta{e) = 
i [5]^(e) — 5]^(e)] describe the lead- molecule coupling. 

For non-local functionals, additional XC contributions 
appear in the lead Fermi functions and give rise to an 



effective bias which is smaller than the physical one [71 
In contrast, the current has Landauer form in the 
ALDA and seems to be equivalent to the static DFT- 
NEGF result. However, the device Green's function 



G£(e) = [el - hu{<7^) - ^l{e) - S^j(e) 



(5) 



is constructed from a Hamiltonian that depends on the 
steady-state density. The latter is not necessarily the 
same as the density obtained from a self-consistent DFT- 
NEGF treatment or even the equilibrium one. To our 
knowledge, a general adiabatic theorem for open systems 
is not available, though adiabaticity has been confirmed 
for model potentials [11]. Such a theorem would prove 
the equivalence of TDLDA and LDA-NEGF approaches 
from Eq. ([5|. Hence, we provide a numerical comparison 
of both approaches for some realistic molecular devices 
in this letter. 

As shown in Ref. ^, the expression for Qa{t) in Eq. ^ 
may be considerably simplified in the adiabatic wide- 
band limit ( AWBL) . It depends finally only on the device 
density matrix ctd, the applied bias V(t) and the WBL 
parameters which characterize the lead self-energies [12] . 
In this way the EOM is closed. At t = 0, the entire 
system is assumed to be in thermal equilibrium and the 
initial density matrix aj:,{t = 0) is easily constructed us- 
ing equilibrium Green's function (EGF) techniques. Sub- 
sequently, bias potentials Va{t) are applied to the con- 
tacts and uniformly shift the lead energy levels. The 
device density matrix is then propagated according 
to Eq. ([T]) , where in each time step the Poisson equation 
is solved to update the device Hamiltonian. 

Results. — We report the calculations of I-V curves for 
two molecular systems, a carbon chain and a benzene- 
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FIG. 3. (a) I-V curves of seven- member ed carbon chain be- 
tween Al nanowires in (001) direction at high bias and low 
bias (inset) regime, (b) TDLDA optical spectrum of isolated 
and hydrogen capped carbon chain (C7H2) (see text). 

diol molecule sandwiched between aluminum leads in the 
(001) direction of bulk Al. The structures are shown 
in Fig. [TJb) and [TJc) respectively. Simulations for an 
Al wire were also performed with similar findings as re- 
ported below. For each of the molecular systems, we in- 
clude explicitly 18 Al atoms of each contact in the device 
region to ensure smooth potentials. The minimal Gaus- 
sian basis set ST0-3G is adopted in the calculations. For 
the propagation, the fourth-order Runge-Kutta method 
is employed with a time step of 0.02 fs to integrate Eq. ([T]) 
in the time domain. The ALDA is adopted for exchange 
and correlation. The bias voltage is turned on expo- 
nentially at t = with a time constant T, according 
to V{t) = Vo{l-e-'/^). 
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FIG. 4. (a) I-V curves of 1,4-benzenediol inbetween Al leads 
at high bias and low bias (inset) regime, (b) TDLDA optical 
spectrum of isolated 1,4-benzenediol (C6O2H6). 

Fig. [2] depicts the transient currents for the two molec- 
ular systems at different values of T with Vb = 3 V. Al- 
though the initial current traces differ, the final steady- 



state current is the same in all cases. This was not 
obvious, as the history of the applied bias voltage is 
included in the formalism, although the memory effect 
of the exchange-correlation functional is absent in the 
ALDA. The establishment of the steady state is ultrafast 
{c^ 10 fs) and occurs in every case. We do not observe per- 
sistent current oscillations as the term in Eq. ([2| 
ensures proper dissipation in the leads. 

TDLDA simulations were then performed for several 
bias values in the linear response, low voltage, regime. 
The asymptotic values of the transient currents (^^lda) 
were extracted in each case and compared to static DFT 
calculations in the LDA using the same WBL approxima- 
tion for the leads in the insets of Fig. |3|a) and Fig.Qa). 
Two kinds of static DFT results are reported. In the first 
case, the device Green's function is determined fully self- 
consistently using the non-equilibrium density (^lda^)- 
In the second case, the equilibrium density is employed 
(/lda) Poisson equation is solved only once, 

which yields a linear potential drop across the molecule. 
For both studied systems all approaches yield identical 
results. This provides numerical evidence for the earlier 
analytical arguments [7 on the equivalence of static and 
dynamic local DFT approaches in linear response. 

In another set of simulations, the covered bias range 
was extended to 5 V. Finite bias simulations allow for the 
coverage of the resonant tunneling regime in which the 
bias window contains one or more of the molecular levels. 
Inspection of the transmission function (not shown here) 
shows that transport occurs predominantly through the 
unoccupied levels in both systems. For benzenediol, the 
conducting molecular orbitals (MO) closest to the Fermi 
energy (Ep) are located 0.92 eV above and 2.58 eV be- 
low Ep^ respectively. The carbon chain exhibits metallic 
behavior with a partially filled lowest unoccupied MO 
derived level [13]. 

The full I-V characteristics of both systems are summa- 
rized in Fig.[3];a) and Fig.[4];a). The TDLDA (Itblba) 
and LDA-NEGF (^lda^) current traces are virtually 
identical over the whole bias range including the region 
of resonant transmission. This level of agreement sup- 
ports not only the validity of the Landauer picture in 
the context of TDLDA (Eq. Q), but also shows that the 
TDLDA and LDA-NEGF densities are in fact identical 
for realistic devices even for strong bias voltage. Another 
reason to study the high bias regime is the possible oc- 
currence of multiple solutions in the self-consistent LDA- 
NEGF treatment. As discussed by Sanchez et al. [9], 
a dynamical approach would single out the most stable 
steady-state solution in such a situation. Attempts to 
find such multiple solutions were unsuccessful for both 
devices. We also performed simulations with increased 
molecule-lead distance in order to cover a larger parame- 
ter space of coupling matrix elements. The transient cur- 
rents show oscillations with slower decay to the steady- 
state in these cases, but again, the TDLDA results match 
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the corresponding LDA-NEGF simulations exactly. 

Comparing now the TDLDA or LDA-NEGF results 
with LDA-EGF values and focussing first on the bias re- 
gion below 2 V, we find significant differences which are 
more pronounced in the carbon chain. The LDA-EGF 
current is initially higher, which can be traced back to a 
movement of the LUMO resonance to higher energies in 
the self-consistent approaches. In fact, TDDFT explicitly 
includes the potential change due to the induced density, 
similar to the random-phase-approximation. In the con- 
text of harmonic perturbations in the optical range, this 
leads to the well known shift of TDDFT excited states 
away from simple KS single-particle energy differences. 
In order to quantify this effect, we computed the optical 
spectrum of the isolated benzenediol and carbon chain 
using a linear response TDDFT method [14. Fig. [3];b) 
andjljb) depict the results using the full response and a 
non-self-consistent treatment employing the ground state 
density, i.e. the KS response. Although in both systems 
the KS and TDDFT spectra differ strongly, the devi- 
ations are higher in the carbon chain with a shift of 
roughly 3 eV compared to a shift of 0.5 eV for benzene- 
diol. This finding is in line with the transport results 
for V < 2 V. For larger values of the bias, the observed 
current traces are more difficult to interpret since the 
transmission profile is altered in a non-linear fashion with 
respect to the equilibrium situation. 

Summary. — In this letter, we performed TDLDA 
transport simulations based on a time propagation of the 
KS density matrix and compared these to conventional 
static DFT results in the Landauer picture. The studied 
systems comprise the ubiquitous benzenediol molecule 
with small equilibrium conductance and an atomic chain 
with metallic conduction properties. In the regime of 
linear response, the analytically predicted equivalence of 
static and time-dependent DFT formalisms for local xc 
potentials could be confirmed with high numerical preci- 
sion. Beyond linear response, TDLDA reproduces the re- 
sult of LDA-NEGF calculations and deviates from static 
LDA only when the equilibrium density is employed in 
the latter. This parallels the energy shifts seen in go- 
ing from the KS response to the full coupled response in 
DFT absorption spectra. Although the transient current 
depends on the specifics of applied bias voltage, a unique 
steady state is reached for the same eventual bias voltage. 

More advanced xc functionals beyond the ALDA are 
expected to lead to an improvement in various ways. For 
example, non-local functionals provide sizable dynami- 
cal corrections to the Landauer formula [7 . Moreover, 
functionals which exhibit a proper derivative disconti- 
nuity refine not only the energetical position of molecu- 
lar levels, but also the molecule-lead coupling [15 . De- 
spite of these short-comings, TDLDA is still useful in the 
qualitative description of transport problems which re- 
quire a dynamical treatment, like the transient behavior 
of molecular devices or photo- switches [16j . 
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